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We show that a single change in the derivation of the linearized semiclassical-initial value representation 
(LSC-IVR or ‘classical Wigner approximation’) results in a classical dynamics which conserves the quantum 
Boltzmann distribution. We rederive the (standard) LSC-IVR approach by writing the (exact) quantum 
time-correlation function in terms of the normal modes of a free ring-polymer (i.e. a discrete imaginary-time 
Feynman path), taking the limit that the number of polymer beads N —)■ oo, such that the lowest normal-mode 
frequencies take their ‘Matsubara’ values. The change we propose is to truncate the quantum Liouvillian, not 
explicitly in powers of at (which gives back the standard LSC-IVR approximation), but in the normal¬ 
mode derivatives corresponding to the lowest Matsubara frequencies. The resulting ‘Matsubara’ dynamics is 
inherently classical (since all terms Olhf) disappear from the Matsubara Liouvillian in the limit N —>■ oo), and 
conserves the quantum Boltzmann distribution because the Matsubara Hamiltonian is symmetric with respect 
to imaginary-time translation. Numerical tests show that the Matsubara approximation to the quantum time- 
correlation function converges with respect to the number of modes, and gives better agreement than LSC- 
IVR with the exact quantum result. Matsubara dynamics is too computationally expensive to be applied to 
complex systems, but its further approximation may lead to practical methods. Copyright (2015) American 
Institute of Physics. This article may be downloaded for personal use only. Any other use requires prior 
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I. INTRODUCTION 


Dynamical properties at thermal equ ilibr ium are of 
central importance to chemical physics.^^^ Sometimes 
these properties can be simulated adequately by entirely 
classical means. But there are plenty of cases, e.g. 
the sp ectr um of liquid water,l^lt^ hydrogen-diffusion on 
metals,!^ and proton/hydride-transfer reactions, ISESl for 
which one needs to evaluate time-correlation functions of 
the form 


|cAB(t) = |Tr 




( 1 ) 


(where Z is the partition function,li^ j3 = I/ZcbT, Tr in¬ 
dicates a complete sum over states, and the other nota¬ 
tion is defined in Sec. II). Such time-correlation functions 
are already approximate, since they employ the quan¬ 
tum Boltzmann distribution (Z in place of the ex¬ 
act quantum-exchange statistics; but this approximation 
is usually adequate (since the thermal wavelength is typi¬ 
cally much smaller than the separations between identical 
particles). What is less well understood is the extent to 
which such functions can be further approximated by re¬ 
placing the exact quantum dynamics by classical dynam¬ 
ics (whilst retaining the quantum Boltzmann statistics). 

The standard way to make this approximation is to 
use the linearized semiclassical-initial value representa¬ 
tion (LSC-IVR, sometimes called the ‘classical Wigner’ 
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approximation) in which the quantum Liouvillian 

is exp ande d as a power series in S^then truncated at hP. 
MiiieiESIIIl and later Shi and Geve!^ showed that this ap¬ 
proximation is equivalent to linearizing the displacement 
between forward and backward Feynman paths in the 
exact quantum time-propagation, which removes the co¬ 
herences, thus making the dynamics classical. The LSC- 
IVR retains the Boltzmann quantum statistics inside a 
Wigner transform,!^ is exact in the zero-time, harmonic 
and high-temperature limits, and has been developed into 
a practical method by several authors.^^SHlSl However, it 
has a serious drawback: the classical dynamics does not 
in general preserve the quantum Boltzmann distribution, 
and thus the quality of the statistics deteriorates over 
time. 

A number of methods have been developed to get 
round this problem, all of which appear to some extent to 
be ad hoc. Some of these methods are obtained by replac¬ 
ing the plain Newtonian dynamics in the LSC-IVR by an 
effective (classical) dynamics which preserves the Boltz¬ 
mann distribution.BSI^ Others, s uch a s the popular cen¬ 
troid molecular dynamics (CMD)^^^ and ring-polymer 
molecular dynamics (RPMD) JHHfelMEll are more heuris¬ 
tic (and still not fully understood) but have the advan¬ 
tage that they can be implemented directly in classi¬ 
cal molecular dynamics codes. An intriguing property 
of CMD and RPMD is that, for some m odel s ystems 
(e.g. the one-dimensional quartic oscillatoi l^^l^^ l), these 
methods give better agreement than LSC-IVR with the 
exact quantum result, even though, like LSC-IVR, they 
completely neglect real-time quantum coherence. 

This last point suggests that the failure of LSC-IVR to 
preserve the quantum Boltzmann distribution may arise, 
not from its neglect of quantum coherence, but from its 
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inclusion of ‘rogue’ components in the classical dynamics. 
The present paper develops a theory that supports this 
speculation. We isolate a core, Boltzmann-conserving, 
classical dynamics, which we call ‘Matsubara dynamics’ 
(for reasons to be made clear). Matsubara dynamics is 
far too expensive to be used as a practical method, but 
is likely to prove useful in understanding methods such 
as CMD and RPMD, and perhaps in developing new ap¬ 
proximate methods. 

The paper is structured as follows. Section II gives key 
background material including the well known ‘Moyal se¬ 
ries’ derivation of the LSC-IVR. Section III re-expresses 
the standard results of Sec. II in terms of ‘ring-polymer’ 
coordinates, involving points along the imaginary-time 
path-integrals that describe the quantum Boltzmann 
statistics. Section IV gives the new results, showing that 
smooth Fourier-transformed combinations of the ring- 
polymer coordinates lead to an inherently classieal dy¬ 
namics which is quantum-Boltzmann-conserving. Section 
V reports numerical tests on one-dimensional models. 
Section VI concludes the article. 


II. BACKGROUND THEORY 

We start by defining the terms and notation to be used 
in classical and quantum Boltzmann time-correlation 
functions (IIA and IIB), and by writing out the standard 
Moyal-series derivation of the LSC-IVR (HC). 


A. Classical correlation functions 

Without loss of generality, we can consider an F- 
dimensional Cartesian system with position coordinates 
q = qi,... ,qp, momenta p, mass m and Hamiltonian 

^(p,q) = + ^(q) (2) 

The thermal time-correlation function between observ¬ 
ables A(p, q), B{p,q) is then 

X A(p,q)H(pt,qt) (3) 

where f dp = dpi ... dpp (and similarly for q), 
and pt = pt (p, q, t) and q* = qt (p, q, t) are the momenta 
and positions after the classical dynamics has evolved for 
a time t. 

Alternatively, we can express B(jpt,qt) as a function 
of the initial phase-space coordinates (p,q): 

H(pt,qt) = R[pt(p,q,t),qt(p,q,t)] = R(p,q,t) (4) 


such that 

X A(p,q)H(p,q,t) 

=_ ^ _ f dn f da 

X A(p,q)e^-*R(p,q,0) (5) 

where the (classical) Liouvillian Cp i^^ 

Tf = ^P • Vq - U(q)Vq • (6) 

with 

dqp ' 

and the arrows indicate the direction in which the deriva¬ 
tive operator is applied (and the backward arrow indi¬ 
cates that the derivative is taken only of U(q) —not of 
any terms that may precede V (q) in any integral). Equa¬ 
tion ^ is less practical than Eq. (|^ (which propagates 
individual trajectories rather than the distribution func¬ 
tion B{p, q, t)) but is better for comparison with the ex¬ 
act quantum expression. 

An essential property of the dynamics is that it pre¬ 
serves the (classical) Boltzmann distribution, which fol¬ 
lows because iL(p, q) is a constant of the motion. As a 
result, we can rearrange Eq. (§ as 

X [e"^^‘A(p, q)] B{p, q, 0) (8) 

showing that CAsit) satisfies 

CAsit) = CBA{—t) ( 9 ) 

which is the detailed balance condition. 

B. Quantum correlation functions 

For clarity of presentation, we will derive the results in 
Secs. Ill and IV for a one-dimensional quantum system 
with Hamiltonian H = T + V, kinetic energy operator 
T = p^/2m, potential energy operator V = V{q), posi¬ 
tion and momentum operators q,p, and mass m. How¬ 
ever, the results we derive in Secs. HI and IV are appli¬ 
cable immediately to systems with any number of dimen¬ 
sions (see Sec. IV.D). 

The simplest form of quantum-Boltzmann time- 
correlation function is that given in Eq. Q, but CAB{t) 
is difficult to relate to the classical time-correlation func¬ 
tion CAsit), because it does not satisfy Eq. and is not 
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in general real. We therefore use the Kubo-transformed 
time-correlation functiorP^ 


CAB{t) = Tr iKpiA) giHtihQ^-iHtih 


( 10 ) 


with 


Kp{A) = - dX (11) 

P Jo 


in powers of h'^. We start by rewriting Eq. (10) as 


CAsit) = dq dA 

J —oo J —oo 

X {q-A/2\Kp{A)\q + A/2) 

x{q + - A/2) (18) 


then insert the momentum identity 


This function gives an equivalent description of the dy¬ 
namics to CABit), to which it is related by a simple 
Fourier-transform formula.^ 

It is easy to show (by noting that and 

commute in Eq. (10)) that CABit) satisfies the detailed 
balance relation 


CABit) = CsAi — t) (12) 

This relation also ensures that CABit) is real (since re¬ 
versing the order of operators in the trace gives Cab it) = 
G*BAi-t))- 

The t = 0 limit of Cab it) can be expressecP^l in terms 
of a classical Boltzmann distribution over an extended 
phase space of ‘ring-polymers’.^^^t^ When A and B are 
functions Aiq) and Biq) of the position operator q, the 
ring-polymer expression is 

X v4(q)B(q)e-^~'^«(P''i) (13) 

where Pn = P/N, fdp = dpi ... dpiv and simi¬ 
larly for f dq, and 

N N 

Biq) = -J^BiqP (14) 

RNiP,^) =TNiP,^)+ UNiq) (15) 

2 ^ 

r»(p,q)=|; + 5^|:fe«-,y (16) 

N 

C/v((q)=^E(ft) (17) 

i=l 

Similar expressions can be obtained when A and B de¬ 
pend on the momentum operator (by inserting position- 
momentum Fourier-transforms). To avoid confusion, we 
emphasise that Eq. © is exact at t = 0, and that we do 
not assume that the ring-polymer Hamiltonian R]\[ip,q) 
generates the dynamics at t > 0. 


1 7°° 

SiA -A') = — dp (19) 

27rh 

to obtain 

-I poo poo 

CABit) =7^ / dq dp 

Zirri J _QQ J—Qo 

X [A:^(i)]w(p, q) [Bit)]wip, q) (20) 
where the Wigner transforms of A and B are given by 

/ OO 

-OO 

x{q-A/2\K^iA)\q + A/2) (21) 


and 

/ OO 

dA 

-OO 

x{q- A/ 2 \C^Gn^ ^/ 2 ). 

( 22 ) 


(and note that we will often suppress the (p, q) depen¬ 
dence of [KpiA)]y^ and [ii(t)jw)- 

We then differentiate Eq. (p0| with respect to t. 


dCABjt) 

dt 



X [KpiA)]y^ 


"^[H.Bit)] 


(23) 


and expand the potential-energy operator in the commu¬ 
tator in powers of A to obtain 


= / dA 

w z —oo 

xiiq-Al2\Bit)\q +AI2) (24) 

with 

mdqdA h ^ X\ dq^ \2 I ^ ^ 

^ A=l,odd ^ ^ ^ 


'j^[H,Bit)] 


C. The LSC-IVR approximation 

1. The Wigner-Moyal series 

To derive the LSC-IVR approximation to CABit), we 
follow ref. 1261 expanding the exact quantum Liouvillian 


Noting that each power of A can be generated by an 
application of —ihd/dp, we then obtain 


dCABjt) _ 1 
dt 2Trh 


dq / dp 


X [Arp(A)]w T[.B(t)]w 


(26) 
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with 


V 1 flhY-^ d^Vjq) 
mdq ^ X\[ 2 J dq>' dp>' ^ ’ 

A=l,odd 


This is the Moyal expansion of the quantum Liouvillian in 
powers of h?. If all terms are included in the series, then 
the application of L generates the exact quantum dynam¬ 
ics (as is easily proved by working backwards through the 
derivation just given). A compact representation of L, 
which will be useful later on is, 


L 

m dq 


V{q)-sm 


ydq 2 dp j ' 


(28) 


(where the paths become infinitessimally short), in the 
harmonic limit (where the are no terms 0{h?) in L), and 
in the high temperature limit (where fluctuations in A 
efficiently dephase) 

Despite these positive features, LSC-IVR suffers from 
the major drawback of not preserving the quantum Boltz¬ 
mann distribution (except in one of the special limits just 
mentioned), since in general 

£[e-/5^]w 0 (33) 


As a result. 




(34) 


where the arrows are defined in the same way as in Eq. ([^ i-®- LSC IVR does not satisfy detailed balance. In 

Secs. III-V we will investigate why this is so. 


2. Approximating the dynamics 


To obtain the LSC-IVR one notes that Eq. (27) can be 
written 


L=C + 0{h‘^) 

where C is the classical Liouvillian 

p d dV d 


C = 


m dq dq dp 


(29) 


(30) 


and then truncates L at The LSC-IVR thus amounts 
to replacing the quantum dynamics by classical dynam¬ 
ics, such that CAB{t) is approximated by 


= / dq dp 


X [Ar/3(A)]w(p, (?) e [B{0)]w{p,q) (31) 


or equivalently 

/ OO nOO 

dq / dp 

-OO j —OO 


X [Ar/3(A)]w(p,q) [i?(0)]w(pt,9t) (32) 

where {pt, qt) are the (classical) position and momentum 
at time t of a trajectory initiated at (p, q) at t = 0. 

Physical insight into the LSC-IVR is obtained by go¬ 
ing back to Eq. (25), and noting that truncating L at 


is equivalent to truncating ( at A. Since A is the 
difference between the origin of a forward path that ter¬ 
minates at z (at time t) and the terminus of a backward 
path that originates at z, it follows that truncating at 
A is equivalent to linearizing the difference between the 
forward and backward Feynman paths at each time-step. 
Hence the neglect of terms 0{h^) is valid if the forward 
and backward paths are very close together, in which 
case there are no coherence effects, and the dynamics be¬ 
comes classical. The LSC-IVR is thus exact at t = 0 


III. RING-POLYMER COORDINATES 

We now recast the standard expressions of Sec. II in 
terms of ring-polymer coordinates. No new approxima¬ 
tions are obtained, but the ring-polymer versions of these 
expressions are needed for use in Sec. IV, where they will 
be used to derive the quantum-Boltzmann-conserving 
‘Matsubara’ dynamics. 


A. Ring-polymer representation of 
Kubo-transformed time-correlation functions 

1. Exact quantum time-correlation function 

Following ref. SSI (see also refs. [T7I and STI) . we define 
the ring-polymer quantum time-correlation function to 
be 




X - Ai_i/2\e-^^^\qi + Ai/2) 

1^1 

X {qi+Ail2\e-^^^/^\zi) 

x{ziW^^/%i-Ai/2) (35) 


where the functions A(q) and B{z) (with z in place of 
q) are defined in Eq. (14) (and we have assumed that 


A and B are functions of position operators to simplify 
the algebra—see Sec. IVD). It is easy to show (by noting 
that A — 1 of all the forward-backward propagators are 
identities, and that the sums in A(q) and B{z) become 
integrals in the limit A —>■ oo) that 


CABit) = lim 

A/—>-oo 


(36) 
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In other words, Eq. (35) in the limit TV —> oo is just 
an alternative way of writing out the standard Kubo- 
transformed time-correlation function CAB{t)- The ad¬ 
vantage of Eq. (35) is that it emphasises the symmetry of 
the entire path-integral expression with respect to cyclic 
permutations of the coordinates qi —>■ gj+i (see Fig. 1); 
this symmetry is otherwise hidden in the conventional 
expression for CAsit) [Eq. (10)]. 


2. Ring-polymer representation of the 
LSC-IVR 


FIG. 1. Schematic diagram showing the structure of the 
(exact) Kubo-transformed quantum time-correlation function 
when represented in ring-polymer coordinates as in Eq. ( |35[ ). 
The red and blue dots represent the coordinates qi and Zi; 
solid lines represent stretches of imaginary time of length divfi; 
arrows represent forward-backward propagations in real time. 


It is straightforward to derive the LSC-IVR approx¬ 
imation from Eq. (35) by generalizing the steps in 


Sec. lie. We insert an identity 


1 

SiAi -A[) = — dpi e^pdA^-^^)/n ( 37 ) 


for each value Z = 1,..., TV, to obtain 


AN]r' 1 


c'AMt) = 


(27rh) 


N 


dq / dp 




_(P,q) B{t) (p,q) 

N L J AT 


(38) 


where 


N 


-PHA 


J N 


(p, q) = / dA A(q) n(9i-i - Ai.,/2\e-^-^\qi + 


(39) 


Z=1 


and 


N 


B(t)J^(p,q)= / dA / dz Biz)]J{qi-Ai/2\e-^^BAzi){zi\e^^^/% + Ai/2)e^P‘^^/^ 


(40) 


1=1 


are generalized Wigner transforms (and we will often sup¬ 
press the dependence on (p,q) in what follows). Note 
that [-jiv and [-j;^ have different forms: [-jAr is a sum of 
products of one-dimensional Wigner transforms, whereas 
[•];jY i® more complicated, with each product coupling 
variables in I and I + 1.1^ Note that since we have speci¬ 
fied that i? is a function of just the position operator (in 
order to simplify the algebra—see Sec. IVD), it follows 
that 


^( 0 ) (p,q) = R(q) 

J N 


(41) 


The next step is to obtain the ring-polymer representa¬ 
tion of the (exact) quantum Liouvillian, which involves 
a straightforward generalization of Eqs. (23)-(27). We 
differentiate C'^(t) with respect to t, obtain a sum of 
N Heisenberg time-derivatives, and expand each mem¬ 
ber in powers of A; to obtain an A-fold generalization of 
Eqs. (24) and (25). On replacing powers of A; by powers 


of —ihd/dpi, we obtain 


dt 


(27r/i) 


N 


dq / dp 




— Ln 

J N 


m 


N 


where 


N 


^ ) h 1 rin, 9 fin, 


1=1 


dqi 2 dpi 


(42) 


(43) 


and the arrow notation is as used in Eq. (H- We can 
write this expression more compactly in terms of t/Ar(q) 
in Eq. 0 as 


Ln = ^P • Vq - 17Ar(q)^ sin^ ^ Vq • 


(44) 
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(since all mixed derivatives of Un{ci) are zero). 

Following Sec. IIC, we then truncate the exact Liou- 
villian at tP such that 


LN=CN + 0{n'^) (45) 


with 


N 


c 


N 


< ^ m rtrii 




m dqi dqi dpi 


(46) 


The ring-polymer version of LSC-IVR thus approximates 
the exact dynamics by the classical dynamics of N inde¬ 
pendent particles, each initiated at a phase-space point 
(pi^qi). The ring-polymer LSC-IVR time-correlation 
function is 


C 


.W[Af] 

AB 


it) = 


1 


(27rh)^ ^ 

1 


dq / dp 


X 


■(27rh)^ ^ 

X e~^^A 


_gCNt 

N 

dq / dp 


m 


N 


N L 


-B(O) (pt,qt) 


N 


(47) 


where 


B{0) (pt,qt) indicates that this Wigner trans- 
L . J Af 

form takes its t = 0 form, but is expressed as a function 
of the momenta and positions (pt,qt) of the N indepen¬ 
dent particles at time t. It is easy show (by noting that 
one can integrate out 1 of the pi) that 


CAsit) = lim 

N—>-(X> 


(48) 


i.e. that the truncation of Tat at gives the standard 
LSC-IVR approximation in the limit TV —> oo (as would 
be expected, since we have approximated the exact quan¬ 
tum Kubo time-correlation function of Eqs. (351 and (36) 
by truncating the quantum Liouvillian at hr). 


B. Normal mode coordinates 


where 


T, 


In 


TV-V2 

■\/2jN sm{2Trln/N) 
\fijN cos(27rZn/iV) 


n = 0 

n = I,...,(V-l)/2 
n = -(TV-l)/2 

(50) 


and similarly for P„ in terms of p;, and in terms of 
A;. The associated normal frequencies take the form 


/SatA 


sm 


/nTT 

\~N 


(51) 


such that the ring-polymer expression for CAniO) 
[Eq. (13)] can be rewritten as 






(52) 


where the normal-mode expression for the ring-polymer 
Hamiltonian (P,Q) is 


^Ar(P, Q) 


(Af-l)/2 

E 


(Ar-l)/2 



+ Un{CI) 


(53) 

and A(Q), B{Q) and Un{Q) are obtained by making the 
substitution 


qi= ^ TinQn (54) 

n=-(Af-l)/2 

into ^(q), B{q) and 17w(q) of Eqs . ([l4l )-p^. Note the 
definition of the sign of a;„ in Eq. ( |5ip , which results in 
somewhat neater expressions later on. Note also that 
i?jv(P, Q) will not be used to generate the dynamics in 
any of the expressions derived below which, like the dy¬ 
namics of Sec. IIIA, will involve TV independent particles 
unconnected by springs. 


1. Definition 

The advantage of ring-polymer coordinates is that we 
can now transform to sets of global coordinates describing 
collective motion of the individual coordinates {pi, g;, A;). 
The choice of global coordinates is not unique. We will 
find it c onven ient to use the normal modes of a free ring- 
polymer, namely the linear combinations of qi that 
diagonalize r/v(p, q) of Eq. (161. These are simply dis¬ 


crete Fourier transforms, which for odd TV (which we will 
assume, to simplify the algebraP^ 


are 


N 


Q„ = ^T,„g,, n = 0,±l,...,±(TV-l)/2 (49) 


1=1 


2. Time-correlation functions 


It is straightforward to convert Eq. (38) into normal 


mode coordinates using the orthogonal transformations 
in Eq. (54), to obtain 


cfiW = 


(27rTi) 


N 


dP / dQ 


-0HA 


_(P,Q) B{t) (P,Q) (55) 

N L J Af 


where 


(7V-l)/2 


dP = 


/ o 

-( 


dP„ 


(56) 


ra=-(Af-l)/2 
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and f dQ is similarly defined. The generalized Wigner 
transforms in Eq. (551 are obtained using Eq. (|54|) to 


substitute (P,Q,D) for (p,q, A) in Eqs. (391 and (|40[), 
and thus contain products of exp(*P„£)„/HJ^in place of 
exTp{ipiAi/h). At t = 0, one obtains 


B{0) (P,Q) = S(Q) 


J N 


(57) 


where B{Q) is obtained by substituting Q for q in i3(q) 
ofEq. (ra. 

The (exact) quantum dynamics is described by 


dC^Akt) 

dt 


(27r/i) 


N 


dP / dQ 




_Ln 

J N 


B{t) 


N 


(58) 


where the Liouvillian is obtained by expressing 
of Eq. (461 in terms of normal modes, which gives 


Ln = —P • Vq — Un{Q)- sin[ - Vq • 
m n \ 2 


(59) 


in which 17Af(Q) is obtained by substituting Q for q in 
C7,v(q) OfEq. (1^. 

As in Sec. IIlA) the LSC-IVR dynamics is obtained by 
truncating Lpf at hP to give 


(7V-r)/2 

Cn = 

n=-{N-l}/2 


Pn_d_ 

m dQr, 


dUNiQ) d 
dQn dPn 


(60) 


after which one obtains in terms of normal 

modes, which gives the (standard) LSC-IVR result in 
the limit V —>■ oo, according to Eq. (481. Hence all we 


have done in Eqs. (55)-(^60| is to re-express the results 


of Sec. IIIA in terms of normal mode coordinates. The 
advantages of doing this will become clear shortly. 


C. Matsubara modes 

We now consider the M lowest frequency ring-polymer 
normal modes in the limit N —> oo, such that M N. 
The frequencies Wn tend to the values 

uin = lim |n|<(M-l)/2 (61) 

Af-foo pn 

which are often referred to as the ‘Matsubara 
frequencies’,!^ and so we will refer to these M modes 
in the limit N ^ oo as the ‘Matsubara modes’. The 
Matsubara modes have the special property that any su¬ 
perposition of them produces a distribution of the coor¬ 
dinates qi which is a smooth and differentiable function 
of imaginary time r, such that 

qi = q{T), T = l = l,...,N (62) 




N ^ oo 


Matsubara 

modes 





FIG. 2. Schematic diagram showing that superpositions of 
Matsubara modes give distributions of path-integral coordi¬ 
nates qi which are smooth, differentiable functions of imagi¬ 
nary time r. Inclusion of non-Matsubara modes gives jagged 
distributions. 


(see Appendix A). Hence distributions made up of super¬ 
positions of the Matsubara modes resemble the sketch in 
Eig. 2. We will often write the Matsubara modes using 
the notation 

Qn= lim n = 0,±l,...,±(M-l)/2 (63) 

AT-foo i/_/V 

(and similarly for P„, Dn). The extra factor of 
ensures that Qn scales as and converges in the limit 
iV —7> oo; e.g. Qo is the centroid (centre of mass) of 
the smooth distribution (?(t). We will refer to the other 
N — M normal modes as the ‘non-Matsubara modes’. In 
general, these modes give rise to jagged (i.e. non-smooth, 
non-differentiable with respect to r) distributions of qi 
(see Fig. 2).!^ 

Matsubara modes have a long historjEzEHIMIlIlin path- 
integral descriptions of equilibrium properties, since they 
give rise to an alternative ring-polymer expression for 
Cab{^)- If we define 

^isko) = ^JdPjdQ A(Q)i?(Q)e-'5«"(P’Q) 

(64) 
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with 


where the ‘Matsubara Liouvillian’ 


Rm{P,Q)=\ E S + 

(65) 


Um(CI)= lim E Tir^VNQr. 

7V_vrvo /V ^ I 


N 


(M-l)/2 


then 


N—>oo N 

/=1 yn=-(M-l)/2 

aM [(M-1)/2]!2 

C'^b( 0)= lim 4“'(0) 

M—^oo 
M^N 


( 66 ) 

(67) 

( 68 ) 


where this limit indicates that M is allowed to tend to 
infinity, subject to the condition that it is always much 
smaller than N, such that the Q remain Matsubara 
modes. In practice, a good approximation to the exact 
result is reached once W(m-i )/2 exceeds the highest fre¬ 
quency in the potential V {q). Equation (64) is less often 


used nowadays to compute static properties, because the 
convergence with respe ct t o M is typically slower than 
the convergence of Eq. (131 with respect to 

However, Eq. (64) tells us something interesting: The 


Boltzmann factor ensures that only smooth distributions 
of (p, q) survive in CAB{t) at t = 0; but at t > 0, 


the force terms in [Eq- (60)] will, in general, mix in 
an increasing proportion of non-smooth, non-Matsubara 
modes, so that the distributions of (p, q) become increas¬ 
ingly jagged as time evolves. The rate at which this mix¬ 
ing occurs depends on the anharmonicity of the potential 
V{q). In the special case that V{q) is harmonic, there is 
no coupling between different normal modes, so the dis¬ 
tributions in (p,q) remain smooth for all time. In other 
words, smooth distributions in (p, q) are found in two of 
the limits (zero-time and harmonic) in which the LSC- 
IVR is known to be exact. 


IV. MATSUBARA DYNAMICS 
A. Definition 


The results of Sec. IIIC suggest that there may be a 
connection between smoothness in imaginary time and 
classical dynamics. We now investigate what happens if 
we constrain an initially smooth function of phase space 
coordinates (p, q) to remain smooth for all (real) times 
t > 0. We take the (exact) quantum Lio uvillian L^v, and 
instead of truncating at hP as in Eq. (601 (which gives the 
LSC-IVR), we retain all powers of ti^take the N —>■ oo 
limit, and split Tjv into 


lim Ln = Cm + lim Lerror(iV, M) (69) 

N—>-oo N—¥oo 


Cm = lim 


Pn d 


(M-l)/2 

y_ 

=-(M-l)/2 ™ 


- C/Ar(Q)^sin| 


(M-l)/2 

E 

n=-(M-l)/2 



(70) 


contains all terms in which the derivatives involve only 
the Matsubara modes, and Terror (iV, M) contains the rest 
of the terms (given in Appendix B). We then discard 
Terror(iV,M), approximating T^r by Cm- We will re¬ 
fer to the (approximate) dynamics generated by Cm as 
‘Matsubara dynamics’. By construction, Matsubara dy¬ 
namics ensures that a distribution of (p, q) which is a 
smooth and differentiable function of r at t = 0 will re¬ 
main so for alH > 0. 

The time-correlation function corresponding to Mat¬ 
subara dynamics is 




N^oo (27rh)^ 


dP / dQ 






J N 


m 


N 


(71) 


We can obtain an explicit form for c'^ab (0 taking the 
same limit as in Eq. (68), allowing M to tend to infinity, 
subject to M <C iV, which gives (see Appendix C) 


/^Mats 

^AB 


it) 


hm cL J (t) 

M^N 


(72) 


where 


= ^j dPj dQ AiQ) 

X e^"*T(Q) 

in which the Matsubara Hamiltonian is 


,-/3[ffAf(P,Q)-ie„(P,Q)] 


(73) 


Hm{P,Q) ^ 7^ + UMiQ) 

2m 


and the phase factor is 


(M-l)/2 

0Af(P,Q)= E Pnl^nQ-r 

n=-(M-l)/2 


(74) 


(75) 


with aM, P and Q defined in Sec. IIIC. Note that, in 
deriving these equations (in Appendix C), we have not 
proved that C'j^^(t) converges with M for t > 0 (only 
that the form of Eqs. (73)-(75) converges with M). We 
test this convergence numerically in Sec. V. 

Thus when the exact dynamics is approximated by 
Matsubara dynamics, the quantum Boltzmann distribu¬ 
tion takes the simple form of a classical Boltzmann dis¬ 
tribution multiplied by a phase factor. At t = 0, one 


















may analytically continue the phase factor (by making 
Pn —>■ Pn — imuJnQ-n) to recover the ring-polymer dis¬ 
tribution in Eq. (64|. However, it is not known whether 
this analytic continuation is valid at t > 0 (except for 
the special case of the harmonic oscillator), and hence 
the most general form of quantum Boltzmann distribu¬ 
tion (in the space of Matsubara modes) is the one given 
in Eq. ([7^. 


B. Matsubara dynamics is classical 


We now rewrite Cm in terms of (P,Q), to make ex¬ 
plicit its dependence on N^ and we also assume that M 


is sufficiently large that Eq. (73) holds, allowing us to 


replace Un{CI)/N by Um{Q)- This gives 

Cm = lim —P • V 5 - 17^m(Q)^ sin( V 5 • ^ 
Af-s-oo m ^ h \2N ^ 


(76) 


In other words, the Moyal series in Matsubara spac^^is 
an expansion in terms of (h/N)"^, rather than Now, 
it is well knowrP^ that the smallness of h cannot in gen¬ 
eral be used to justify truncating the (standard LSC- 
IVR) Moyal series of Eq. (27) at since at least one of 
the Wign er transforms in the time-correlation function 
[Eq. ( [20| ] contains derivatives that scale as How¬ 
ever, it is easy to show that the derivatives of all terms 
in the integral in Eq. (73) scale as . As a result, it 


follows that all derivatives higher than first order in Cm 
vanish in the limit N ^ 00 , with the result that 


C 


M 


= E 

n=-(M-l)/2 


Pn d _ dllMiO.) d 
dQn dQn dPn 


(77) 


In other words, Matsubara dynamics is classical. 

This is a surprising result, which needs to be inter¬ 
preted with caution. It does not mean that the depen¬ 
dence of B{Q) on the Matsubara modes evolves classi¬ 
cally in the exact quantum dynamics, since the exact 
Liouvillian Lj^ contains derivative terms that couple the 
Matsubara modes with the non-Matsubara modes (for 
which the higher-order derivatives cannot be neglected): 
it means that the dynamics of the Matsubara modes be¬ 
comes classical when they are decoupled from the non- 
Matsubara modes. 

One way to understand the origin of the h/N in 
Eq. (76) is to note that the Fourier transform between 
Pn and Dn (in the Wigner transforms of Eqs. (39) and 


(40)) is exp{iNPnDn/li). Hence the effective Planck’s 
constant associated with motion in the Matsubara coor¬ 
dinates tends to zero in the limit N ^ 00 . Note that 
the dependence of the Boltzmann distribution on the 
non- Mat subara modes is more complicated than that of 
Eq. (73), and contains powers of {h/N)~^ which cancel 


out the powers of (fi/N) in Ljv (which must obviously 


happen, since we know that the exact dynamics is not in 
general classical). 

Matsubara dynamics thus has many features in com¬ 
mon with LSC-IVR: it is exact in the t = 0 limit (when 
all distributions of (p,q) are smooth superpositions of 
Matsubara modes), in the harmonic limit (where the dy¬ 
namics of the Matsubara modes is decoupled from that 
of the non-Matsubara modes), and in the classical limit 
(since setting M = 0 in Eq. (73) gives the classical time- 
correlation function); and it neglects all terms 0{h^) 
in the (exact) quantum Liouvillian. However, Matsub¬ 
ara dynamics differs from LSC-IVR in that it also ne¬ 
glects the terms 0{1tP) that contain derivatives in the 
non-Matsubara modes. One can thus regard Matsubara 
dynamics as a filtered version of LSC-IVR, in which the 
parts of the dynamics that cause the smooth distribu¬ 
tions of (p,q) to become jagged have been removed.^ 


C. Conservation of the quantum Boltzmann 
distribution 


Confining the dynamics to the space of Matsubara 
modes has a major effect on the symmetry of the Hamil¬ 
tonian. The LSC-IVR Hamiltonian i7jv(Pj Q) is simply 
the classical Hamiltonian of N independent particles, and 
is thus symmetric with respect to any permutation of the 
phase space coordinates [e.g. {pi,qi) -n- (^ 3 , 93 )]. On re¬ 
stricting the dynamics to the Matsubara modes, most 
of these symmetries are lost (since individual permuta¬ 
tions would destroy the smoothness of the distributions 
of (p,q))- However, one operation which is retainecP^ 
is symmetry with respect to cyclic permutation of the 
coordinates, which, on restricting the dynamics to Mat¬ 
subara space, becomes a continuous, differentiable sym¬ 
metry, namely invariance with respect to translation in 
imaginary time: 


dgM(P,Q) 

dr 


= 0 


(78) 


(see Appendix A). It thus follows from Noether’s 
theorem,that 


(iAM(P,Q) ^ d^ 

dr dt 


(M-l)/2 

E ■ 

I n=-(M-l)/2 


dQn 

dr 


= 0 (79) 


where Am(P, Q) is the Matsubara Lagrangian. In other 
words, in Matsubara dynamics, there exists a constant 
of the motion (in addition to the total energy) which is 
given by the term in brackets above. 

In Appendix A, it is shown that the phase 6m{P, Q) 
in the quantum Boltzmann distribution [Eqs. (73)- 
can be written 


(M-l)/2 

9m(P,Q) = - Y. 

ra=-(M-l)/2 


dQn 

dr 


(80) 
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and is thus the consterh of the motion associated with the 
invariance of Hm(P,Q) to imaginary time-translation. 
Since Hm (P,Q) is of course also a constant of the mo¬ 
tion, it follows that Matsubara dynamics conserves the 
quantum Boltzmann distribution. 

As a result, Matsubara dynamics satisfies the detailed 
balance relation 


( 81 ) 

and gives expectation values 

=(B)''^'(0) (82) 


which are independent of time (and equal to the exact 
quantum result in the limit M —)■ oo; see Eq. (68)). Note 
that the step between the second and third lines follows 
from analytic continuation (P„ —> — imUnQ-n)- 

We thus have the surprising result that a purely clas¬ 
sical dynamics (Matsubara dynamics) which uses the 
smoothed Hamiltonian that arises naturally when the 
space is restricted to Matsubara modes, conserves the 
quantum Boltzmann distribution. At first sight this may 
appear counter-intuitive. For example, it is clear that the 
classical dynamics will not respect zero-point energy con¬ 
straints, nor will it be capable of tunnelling. However, it 
is the phase 6m{P, Q) which converts what would other¬ 
wise be a classical Boltzmann distribution in an extended 
phase-space into a quantum Boltzmann distribution, and 
the phase is conserved. 


D. Generalizations 

The derivations above can easily be generalized to sys¬ 
tems with any number of dimensions. For a system whose 
classical Hamiltonian resembles Eq. ([^ , there are F x M 
Matsubara modes, one set of M modes in each dimension. 
All the steps in Secs. HI and IV.A-C are then the same, 
except that, with P dimensions instead of one, there is 
now a sum of P phase terms, each resembling 0 m(P, Q)- 
Noether’s theorem shows that the sum of these terms and 
hence the quantum Boltzmann distribution is conserved. 

We emphasise that the derivations above were carried 
out for operators A and B in CAB{t) which are general 
functions of the coordinate operators q. Matsubara dy¬ 
namics is therefore not limited to correlation functions 



t (a.u.) 

FIG. 3. Convergence with respect to number of modes M 
of the Matsubara position auto-correlation function 
calculated for the quartic potential of Eq. ( |83[ ), at a reciprocal 
temperature of /3 = 2 a.u. The red lines correspond to M = 1 
(dots), 3 (chains), 5 (dashes) and 7 (solid). The solid black 
line is the exact quantum result. 



t (a.u.) 

FIG. 4. Evolution of the phase 6m{P, Q) along a single clas¬ 
sical trajectory on the quartic potential, with M = 5, and 
N = 5 (dots), 9 (dashes) and oo (solid line). The latter 
corresponds to Matsubara dynamics in which the phase is 
conserved. 


involving linear operators of position. The derivations 
can also be repeated, with minor modifications in the al¬ 
gebra, for the case that A and B are general functions of 
the momentum operator (which results in functions of P 
appearing in the generalised Wigner transforms). 


V. NUMERICAL TESTS OF THE 
EFFICACY OF MATSUBARA DYNAMICS 

So far we have made no attempt to justify the use of 
Matsubara dynamics, beyond pointing out that it is exact 
in all the limits in which LSC-IVR is exact, but that, un¬ 
like LSC-IVR, it also conserves the quantum Boltzmann 
distribution. Here we investigate whether Matsubara dy¬ 
namics converges with respect to the number of modes 
M, and make numerical comparisons with the LSC-IVR, 
CMD and RPMD methods. 
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FIG. 5. Time-dependence of the thermal expectation value 
for the quartic potential at /3 = 2, computed us¬ 
ing LSC-IVR (blue), and Matsubara dynamics (red: M = 1 
(dots), 3 (dashes), 5 (solid)), and compared with the exact 
quantum result (black). 


FIG. 6. Comparisons of position-autocorrelation functions 
computed using different levels of theory, for (a) the quartic 
potential and (b) the weakly anharmonic potential of Eq. (841. 
The Matsubara results were obtained using M = 7 (quartic) 
and M = 5 (weakly anharmonic). 


The presence of the phase 0 m(P, Q) in the Boltzmann 
distribution [Eq. (73l] means that Matsubara dynamics 
suffers from the sign problem, and thus cannot be used 
as a practical method. Howe ver, we were able to evaluate 
(i.e. of Eq. (73) with A = q, B = q) for 

some one-dimensio nal m odel systems. For consistency 
with previous work,l22lS31 we considered the quartic po¬ 
tential 


V{q) = \q^ 


(83) 


and the weakly anharmonic potential 

, 84 ) 

where atomic units are used with m = 1. Calculations 
using potentials with intermediate levels of anharmonic- 
ity were found to give similar results (and are not shown 
here). 

Figure 3 shows for the quartic potential, at an 

inverse temperature of /3 = 2 a.u., for various values of 
M. These results were obtained by propagating classi¬ 
cal trajectories using the Matsubara potential UmIQ) to 
generate the forces, subject to the Anderson thermostalP 
(according to which each was reassigned to a value 
drawn at random from the classicM Bojtzmann distribu¬ 
tion every 2 atomic time units); Um{Q) was computed 


by taking the TV —> oo limit analytically, as described in 
the supplemental material.l^ A total of 10^^ Monte Carlo 
points was found necessary to converge (t). Extend¬ 
ing these calculations beyond M = 7 was prohibitively 
expensive, and the final few M were particularly diffi¬ 
cult to converge (since 0 m(P,Q) becomes increasingly 
oscillatory as increases). Nevertheless, the results in 
Fig. 3 are sufficient to show that C'q^^(t) converges with 
respect to M, although the convergence appears to be¬ 
come slower as t increases. For the weakly anharmonic 
potential, convergence to within graphical accuracy was 
obtained using M = 5 for /3 = 2 a.u. 

We also conhrmed numerically that Matsubara dynam¬ 
ics conserves the quantum Boltzmann distribution. Fig¬ 
ure 4 shows the phase 0 m(P5 Q) as a function of time 
along a Matsubara trajectory. When a coarse number of 
polymer beads {N = 5) is used, such that the M lowest- 
frequency modes are a poor approximation to the Mat¬ 
subara modes, the phase is not conserved; however, as 
N is increased, the variation of the phase along the tra¬ 
jectory flattens, becoming completely time-independent 
in the limit N —> oo. Figure 5 plots the expectation 


value which is found to be time-independent 

as expected from Eq. (82). 

Figure 6 compares the Matsubara correlation functions 
(t) for both potentials with exact quantum, LSC- 
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IVR, CMD and RPMD results. The quartic potential at 
j3 = 2 (panel a) is a severe test for which any method 
that neglects real-time coherence fails after a single re¬ 
currence. Nevertheless, we see that Matsubara dynamics 
gives a much better treatment than LSC-IVR, reproduc¬ 
ing almost perfectly the first recurrence at 6 a.u., and 
damping to zero more slowly.^S The Matsubara result 
is also better than both the CMD and RPMD results. 
The same trends are found for the weakly anharmonic 
potential (Fig. 6, panel b), and were also found for the 
potentials with intermediate anharmonicity (not shown). 


VI. CONCLUSIONS 

We have found that a single change in the derivation 
of LSC-IVR dynamics gives rise to a classical dynam¬ 
ics (‘Matsubara dynamics’) which preserves the quantum 
Boltzmann distribution. This change involves no explicit 
truncation in powers of but instead a decoupling of 
a subspace of ring-polymer normal modes (the Matsub¬ 
ara modes) from the other modes. The dynamics in this 
restricted space is found to be purely classical and to en¬ 
sure that smooth distributions of phase-space points (as 
a function of imaginary time), which are present in the 
Boltzmann distribution at time t = 0, remain smooth at 
all later times. The LSC-IVR dynamics, by contrast, in¬ 
cludes all the modes, which has the effect of breaking up 
these smooth distributions, and thus failing to preserve 
the quantum Boltzmann distribution. Numerical tests 
show that Matsubara dynamics gives consistently better 
agreement than LSC-IVR with the exact quantum time- 
correlation functions. 

These results suggest that Matsubara dynamics is 
a better way than LSC-IVR, at least in principle, to 
account for the classical mechanics in quantum time- 
correlation functions. We suspect that Matsubara 
dynamics may be equivalent to expanding the time- 
dependence of the quantum time-correlation function in 
powers of and truncating it at this is in con¬ 
trast to LSC-IVR, in which one truncates the quantum 
LiouvilliarP^ at hP. However, further work will be needed 
to prove or disprove this conjecture. 

Matsubara dynamics is far too expensive to be useful 
as a practical method. However, it is probably a good 
starting point from which to make further approxima¬ 
tions in order to develop such methods. The numerical 
tests reported here show that Matsubara dynamics gives 
consistently better results than both CMD and RPMD, 
suggesting that these popular methods may be approxi¬ 
mations to Matsubara dynamics. 
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Appendix A: Differentiability with respect to 
imaginary time 


A distribution of ring-polymer coordinates qi, I = 
1,... ,iV, can be written as a smooth and differentiable 
function of the imaginary time t (0 < r < /3h) if the 
limit 


dq{T) 

dr 


gz+i - qi-i 

lim -—r-, 

n—¥oo 2l3j\ih 


T = Ih^N 


(Al) 


exists, i.e. if 


lim qi+i - q/_i ~ N ^ (A2) 

N—¥CC> 


For a distribution formed by superposing only the Mat¬ 
subara modes, we can use trigonometric identities and 
the definitions in Sec. HI to write 


qi+i - qi-i 

(M-l)/2 I- 


= 2 V 2 


X sm 


E 

n—l 


COS 


/ 2Tml 


Qi 


- sin ( 


2'Knl 


Q-r 


( 27rn\ 

[~w) 


(A3) 


Since n N, the sine function on the right ensures that 
Eq. (A2| is satisfied; also, repetition of this procedure 


shows that higher-order differences of order A scale as 
N~^. Hence a distribution qi formed from a superposi¬ 
tion of Matsubara modes is a smooth and differentiable 
function of t. The same is true for distributions in pi and 

A;. 

To prove that the Matsubara Hamiltonian is invariant 
under imaginary-time translation [Eq. (78)], we first dif¬ 
ferentiate the Matsubara potential Um{Q) with respect 
to T, which gives 


dUMiQ) _ Vi Um( cj) - UM(q) 

dr 00 


(A4) 


where 


N 


' N (M-l)/2 \ 

Um{<1) = E^IE E TlnTmnqm (AS) 

i = l ym=l „=_(M-1)/2 ) 

and Vi represents a cyclic permutation of the coordinates 
q-m qra+h SUCh that 

^ N / N (M-l)/2 

'Pi Um{(i) = E^ E E 'Pln'P(m—l) ndm 

1=1 \m=l n=-{M-l)l2 


TJHH, MJW and SCA acknowledge funding from the 
UK Science and Engineering Research Council. AM ac¬ 
knowledges the European Lifelong Learning Programme 
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We then rearrange the sum over n in Eq. (A6) into 


(M-l)/2 

210^(771-1)0+ [21n+'(m-l) 71 + 2/_nT(m-l)-n] 

n—1 

(A7) 


and use trigonometric identities to show that 

2l7iE(777— 1 ) n + E/ —nT^irn—l) —n 


— T(i + i) nTmn + T( 


(/ + 1) n^mn “t“ (/ + 1) —n^m —n 


(A8) 


Re-ordering the sum over I and using the property that 
Tio = TV-1/2 gives 


Vi c/M(q) = c/M(q) 


which proves that 


dUMiQ) 

dr 


= 0 


(A9) 


(AlO) 


The same line of argument can be applied to the kinetic 
energy 'P'^f2m, thus proving Eq. (78). 

To obtain the derivative of Qn with respect to r 
(needed to prove Eq. (801), we write 

dQ Y] 


1 


N 


= lim Ti, 

dr Af-7oo y/n ^ 

N 


qi +1 - qi-i 
' 2PNh 


= lim ^ [T(l-l)n T{l+l)n](ll 

N^oo y(/V ^ 2/3n1^ 

and use trigonometric identities to obtain 
Ti+in - Ti^in = 2Ti_nSin{2mr/N) 
Since n N, it follows that 
dQji 


dr 


— ^nQ—r. 


(All) 


(A12) 


(A13) 


where ojn is the Matsubara frequency defined in Eq. (61). 


Appendix B: Error term for Matsubara 
Liouvillian 


The error term Terror (A^,+7) of Eq. (69) is the differ¬ 
ence Ln — Cm between the exact quantum Liouvillian 
and the Matsubara Liouvillian. Using Eqs. (59) and (70) 
and the trigonometric identity 


sin(a + b) — sin a = 2 sin [ - 


cos a -I- - 


(Bl) 


we can write 


Terror (TV, M) = 


P-7 


d 


(W-l)/2 

77 = (Si)/ 2 ^^ 


4 X 

2 


m dQr, 



(B2) 


with 


h _^ 

^ n=(M+l)/2 '^'^77 dPn 


(B3) 


and 



(M-l)/2 

E 

n=-(M-l)/2 


dQn dPn 


(B4) 


Appendix C: Derivation of Matsubara 
time-correlation function 


To obtain the expression for in Eq. (73), 


we note that B (P,Q,/) is independent of the non- 
Matsubara P modes (since, by construction, these modes 
are not involved in the Matsubara dynamics) which can 
therefore be integrated out, giving a product of Dirac 
delta-functions in the non-Matsubara D modes.l^_As a 


result, the Wigner transform 
duces to 


-PHa 


_ in Eq. (71) re- 
J AT 


-PHa 


J N 
N-M 


(Pm, Q) 


(M-l)/2 

= (27rT)i''-“A(Q) / (TDm H 

■/ 77=-(M-1)/2 

N 

X n(i?r-i(Q>DM)|e-^-'^|7?+(Q,DM)) (Cl) 
1=1 


where Pm and Dm include only the Matsubara modes 
(and Q includes all N modes), and 


(Af-l)/2 (M-l)/2 

^f(Q,DM)= y]] TinQnP^ y]] TinDn/2 

„=-(JV-l)/2 n=-(M-l)/2 

(C2) 

(where the dependence of rjf on (Q, Dm) will be sup¬ 
pressed in what follows). Expressing the bra-ket in ring- 
polymer form, and using trigonometric identities, we ob¬ 
tain 


,-PHa _(Pm,Q) 
J N 


N/2 


A(Q) / dDM 


(M-l)/2 

X giPnDri/h 


7i=-(M-1)/2 


J Pn 

xexp - — 


N 


) + V{v^) 




(C3) 
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where 

/m(Q,Dm) 

4 


(M-l)/2 

E 


(PNflY 

{N-l)/2 

+ E ^-n) 

n=(M+l)/2 


Qn sin 


riTT D 


N 


-n nn 

cos — 
2 N 


(C4) 


On taking the limit N —> cx), and converting D m to 
D, we find that the Gaussians involving Dm in Eq. ( |C3| 
have the form 


exp {-mDlN'^ 


(C5) 


i.e. each Gaussian in D becomes a Dirac delta-function 
in the limit N — » oo. This allows us to replace the third 
line in Eq. (C3) by 


(C6) 



N 

/ i.N-l)/2 \‘ 

exp 


TlnQn j 



\n=-{N-l)/2 j _ 


and to integrate out the D, giving 




_(Pm, Q) 


N 


V Hn ) 


dl(Q) 


X e 


X exp 


X exp 


(M-l)/2 

—/3ivPM/2m I I 2iP-„,Q-ntan(mr/N)/h 


n - 

ri=-(M-l)/2 


N 


.V I (JV-l)/2 

-TT.y{ T. r,„< 3 , 

1=1 \n=-{N-l)/2 


Pnip 


(N-l)/2 


E (Qn + Q- 


2 

n/^n 


n=(M+l)/2 


(C7) 


We then substitute this expression into the integral of 
Eq. ( [7l| ) (with J dP replaced by / c?Pm)) and take the 
limit iW —>■ oo (subject to M ^ N), which allows us to 
integrate out the non-Matsubara modes in Q. Use of the 
formulcP 


N-1 


sin (mr/N) = N/2 


N-l 


(C8) 


then gives Eqs. (73)-(751. 
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